library(tidyverse)
library(RColorBrewer)
library(readbitmap)
library(rjson)
library(patchwork)
geom_spatial <- function(mapping = NULL,
                         data = NULL,
                         stat = "identity",
                         position = "identity",
                         na.rm = FALSE,
                         show.legend = NA,
                         inherit.aes = FALSE,
                         ...) {
  
  GeomCustom <- ggproto(
    "GeomCustom",
    Geom,
    setup_data = function(self, data, params) {
      data <- ggproto_parent(Geom, self)$setup_data(data, params)
      data
    },
    
    draw_group = function(data, panel_scales, coord) {
      vp <- grid::viewport(x=data$x, y=data$y)
      g <- grid::editGrob(data$grob[[1]], vp=vp)
      ggplot2:::ggname("geom_spatial", g)
    },
    
    required_aes = c("grob","x","y")
    
  )
  
  layer(
    geom = GeomCustom,
    mapping = mapping,
    data = data,
    stat = stat,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}

Introduction

This example shows one how to estimate immunofluorescence intensity at each spot. It is for educational purposes only and is not supported by 10x genomics.

We will be using a public data set from: https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Adult_Mouse_Brain_Coronal_Section_1

Note: all paths are relative to my home directory and will not work on your system. Adjust them as needed.

The image channel setup. This is different from what is on the support site. The support documentation is correct for the loupe file

  • [,,1] = FITC NeuN
  • [,,2] = DAPI DAPI
  • [,,3] = TRITC Fiducial Frame

Get the data

tmpdir <- tempdir()

url <- 'https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain_Coronal_Section_1/V1_Adult_Mouse_Brain_Coronal_Section_1_spatial.tar.gz'
file <- basename(url)
download.file(url, file)

untar(file, exdir = tmpdir)
list.files(tmpdir, recursive = TRUE)
## [1] "spatial/aligned_fiducials.jpg"     "spatial/detected_tissue_image.jpg"
## [3] "spatial/scalefactors_json.json"    "spatial/tissue_hires_image.png"   
## [5] "spatial/tissue_lowres_image.png"   "spatial/tissue_positions_list.csv"
images_cl <- read.bitmap("spatial/tissue_lowres_image.png")
images_chnl_1 <- images_cl[,,1]
images_chnl_2 <- images_cl[,,2]
images_chnl_3 <- images_cl[,,3]

Take a quick look at the image channels

image(images_chnl_1)

image(images_chnl_2)

image(images_chnl_3)

Get the height and width of the images

height <-  data.frame(height = nrow(images_cl))
width <- data.frame(width = ncol(images_cl))

Make the images into grobs for easy plotting in ggplot2

grobs <- list(grid::rasterGrob(images_cl, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_1 <- list(grid::rasterGrob(images_chnl_1, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_2 <- list(grid::rasterGrob(images_chnl_2, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_3 <- list(grid::rasterGrob(images_chnl_3, width=unit(1,"npc"), height=unit(1,"npc")))


images_tibble <- tibble(grob=grobs, 
                        grob_chnl_1 = grobs_chnl_1, 
                        grob_chnl_2 = grobs_chnl_2,
                        grob_chnl_3 = grobs_chnl_3)
images_tibble$height <- height$height
images_tibble$width <- width$width
images_tibble
## # A tibble: 1 x 6
##   grob       grob_chnl_1 grob_chnl_2 grob_chnl_3 height width
##   <list>     <list>      <list>      <list>       <int> <int>
## 1 <rastrgrb> <rastrgrb>  <rastrgrb>  <rastrgrb>     600   600

Get the scale factors so that we can adjust the spot pixel positions for the lowres image

scales <- fromJSON(file = "spatial/scalefactors_json.json")

Make a merged dataframe with all the information we need.

bcs <- read.csv("spatial/tissue_positions_list.csv",
                 col.names=c("barcode","tissue","row","col","imagerow","imagecol"), header = F) 
bcs$imagerow_scaled <- bcs$imagerow * scales$tissue_lowres_scalef    # scale tissue coordinates for lowres image
bcs$imagecol_scaled <- bcs$imagecol * scales$tissue_lowres_scalef
bcs$imagerow_scaled_round <- round(bcs$imagerow * scales$tissue_lowres_scalef) # Rounded scales
bcs$imagecol_scaled_round <- round(bcs$imagecol * scales$tissue_lowres_scalef)
bcs$tissue <- as.factor(bcs$tissue)
bcs$height <- height$height
bcs$width <- width$width

Get intensities for each pixel

Function to get intensities

This is taking into account the spot diameter and pulling the average of pixels in every direction Might be better if using the full resolution image but seems to perform well on the lower resolution images.

get_intensity <- function(bcs, image, scales) {

radius <- round((scales$spot_diameter_fullres * scales$tissue_lowres_scalef)/2)
intensity <- list()

for (i in 1:nrow(bcs)) {
  temp <- expand.grid(x = (bcs$imagerow_scaled_round[i] - radius):(bcs$imagerow_scaled_round[i] + radius), 
                      y = (bcs$imagecol_scaled_round[i] - radius):(bcs$imagecol_scaled_round[i] + radius))
  
  intensity[[i]] <- mean(image[temp$x, temp$y])
  
}
return(intensity)
}

Calculate intensities for each channel and each spot

intensity_1 <- get_intensity(bcs, images_chnl_1, scales)
intensity_2 <- get_intensity(bcs, images_chnl_2, scales)
intensity_3 <- get_intensity(bcs, images_chnl_3, scales)

Add the intensity info to bcs

bcs$intensity_1 <- unlist(intensity_1)
bcs$intensity_2 <- unlist(intensity_2)
bcs$intensity_3 <- unlist(intensity_3)

head(bcs)
##              barcode tissue row col imagerow imagecol imagerow_scaled
## 1 ACGCCTGACACGCGCT-1      0   0   0     2238     3201        55.39604
## 2 TACCGATCCAACACTT-1      0   1   1     2478     3337        61.33663
## 3 ATTAAAGCGGACGAGC-1      0   0   2     2240     3476        55.44554
## 4 GATAAGGGACGATTAG-1      0   1   3     2480     3612        61.38614
## 5 GTGCAAATCACCAATA-1      0   0   4     2242     3751        55.49505
## 6 TGTTGGCTGGCGGAAG-1      0   1   5     2482     3887        61.43564
##   imagecol_scaled imagerow_scaled_round imagecol_scaled_round height width
## 1        79.23267                    55                    79    600   600
## 2        82.59901                    61                    83    600   600
## 3        86.03960                    55                    86    600   600
## 4        89.40594                    61                    89    600   600
## 5        92.84653                    55                    93    600   600
## 6        96.21287                    61                    96    600   600
##   intensity_1 intensity_2 intensity_3
## 1  0.02352941  0.01176471  0.04329412
## 2  0.02352941  0.01176471  0.04313725
## 3  0.02352941  0.01176471  0.04454902
## 4  0.02352941  0.01176471  0.04313725
## 5  0.02352941  0.01176471  0.04313725
## 6  0.02352941  0.01176471  0.04313725

Plots

Define our color palette for plotting.

myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))

I like this color scale but it’s easily changed with code like this

scale_fill_gradient(low = "black", high = "red")+

instead of

scale_fill_gradientn(colours = myPalette(100))+

Plot the data with and without spot intensities

plot1 <- bcs %>% 
  filter(tissue =="1") %>% 
      ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_1)) +
                geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_1), x=0.5, y=0.5)+
                geom_point(shape = 21, colour = "black", size = 2, stroke = 0.1)+
                coord_cartesian(expand=FALSE)+
                scale_fill_gradientn(colours = myPalette(100))+
                xlim(0,max(bcs %>% 
                            dplyr::select(width)))+
                ylim(max(bcs %>% 
                            dplyr::select(height)),0)+
                xlab("") +
                ylab("") +
                labs(fill = "NeuN Intensity")+
                ggtitle("NeuN Spot Intensity") +
                theme_set(theme_bw(base_size = 10))+
                theme(panel.grid.major = element_blank(), 
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(), 
                        axis.line = element_line(colour = "black"),
                        axis.text = element_blank())
plot2 <- bcs %>% 
  filter(tissue =="1") %>% 
      ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_1)) +
                geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_1), x=0.5, y=0.5)+
                coord_cartesian(expand=FALSE)+
                xlim(0,max(bcs %>% 
                            dplyr::select(width)))+
                ylim(max(bcs %>% 
                            dplyr::select(height)),0)+
                xlab("") +
                ylab("") +
                ggtitle("NeuN Image Only") +
                theme_set(theme_bw(base_size = 10))+
                theme(panel.grid.major = element_blank(), 
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(), 
                        axis.line = element_line(colour = "black"),
                        axis.text = element_blank())

plot1 + plot2

plot1 <- bcs %>% 
  filter(tissue =="1") %>% 
      ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_2)) +
                geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_2), x=0.5, y=0.5)+
                geom_point(shape = 21, colour = "black", size = 2, stroke = 0.1)+
                coord_cartesian(expand=FALSE)+
                scale_fill_gradientn(colours = myPalette(100))+
                xlim(0,max(bcs %>% 
                            dplyr::select(width)))+
                ylim(max(bcs %>% 
                            dplyr::select(height)),0)+
                xlab("") +
                ylab("") +
                labs(fill = "DAPI Intensity")+
                ggtitle("DAPI Spot Intensity") +
                theme_set(theme_bw(base_size = 10))+
                theme(panel.grid.major = element_blank(), 
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(), 
                        axis.line = element_line(colour = "black"),
                        axis.text = element_blank())
plot2 <- bcs %>% 
  filter(tissue =="1") %>% 
      ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_2)) +
                geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_2), x=0.5, y=0.5)+
                coord_cartesian(expand=FALSE)+
                xlim(0,max(bcs %>% 
                            dplyr::select(width)))+
                ylim(max(bcs %>% 
                            dplyr::select(height)),0)+
                xlab("") +
                ylab("") +
                ggtitle("DAPI Image Only") +
                theme_set(theme_bw(base_size = 10))+
                theme(panel.grid.major = element_blank(), 
                        panel.grid.minor = element_blank(),
                        panel.background = element_blank(), 
                        axis.line = element_line(colour = "black"),
                        axis.text = element_blank())

plot1 + plot2